import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import os

def RL_calculation(list,d,mode):
    # 计算非磁性材料没有问题,计算磁性材料没有问题
    # d为厚度，单位mm
    c = 299792458
    if mode == 1:
        a = (list[:,3] - list[:,4]*1j) / (list[:,1] - list[:,2]*1j)
        b = (list[:,3] - list[:,4]*1j) * (list[:,1] - list[:,2]*1j)
        z = np.sqrt(a) * np.tanh(1j*2*d*1e-3*np.pi*list[:,0]*np.sqrt(b)/c)
        RL = 20*np.log10(abs((z-1)/(z+1)))
        Zin_Z0 = abs(z)
    else:
        a = 1 / (list[:,1] - list[:,2]*1j)
        b = 1 * (list[:,1] - list[:,2]*1j)
        z = np.sqrt(a) * np.tanh(1j*2*d*1e-3*np.pi*list[:,0]*np.sqrt(b)/c)
        RL = 20*np.log10(abs((z-1)/(z+1)))
        Zin_Z0 = abs(z)
    return RL, Zin_Z0

def attenuation_calculation(list, mod):
    # 计算磁性材料没有问题，计算非磁性材料没有问题
    c = 299792458
    if mod == 1:
        para_1 = list[:,4]*list[:,2] - list[:,3]*list[:,1]   # u''e''-u'e'
        para_2 = list[:,3]*list[:,2] + list[:,4]*list[:,1]   # u'e''+u''e
        sqrt_1 = np.sqrt(para_1**2 + para_2**2)
        sqrt_2 = np.sqrt(para_1 + sqrt_1)
        attenuation = np.sqrt(2)*np.pi*list[:,0]*sqrt_2/c
    else:
        para_1 = 0*list[:,2] - 1*list[:,1]   # u''e''-u'e'
        para_2 = 1*list[:,2] + 0*list[:,1]   # u'e''+u''e
        sqrt_1 = np.sqrt(para_1**2 + para_2**2)
        sqrt_2 = np.sqrt(para_1 + sqrt_1)
        attenuation = np.sqrt(2)*np.pi*list[:,0]*sqrt_2/c
    return attenuation

def data_conversion(list):
    """
    将一个电磁参数数据变为标准的5列数据,依次是频率，介电常数实部和虚部，磁导率实部和虚部
    将频率变为Hz
    :param list: 电磁参数数据，列数应为3或5
    :return: 转换后的数据，包括频率、介电常数实部和虚部、磁导率实部和虚部，共5列
    """
    data_row, data_col = list.shape
    # 频率, Hz
    if list[0, 0] > 1e9:
        frequency = list[:, 0]
    else:
        frequency = list[:, 0] * 1e9

    # 介电常数实部和虚部
    e_1 = list[:, 1]
    e_2 = list[:, 2]

    # 磁导率实部和虚部
    if data_col == 5:
        u_1 = list[:, 3]
        u_2 = list[:, 4]
    else:
        u_1 = np.ones(data_row)
        u_2 = np.zeros(data_row)
    # 结果
    data_conversion = np.column_stack((frequency, e_1, e_2, u_1, u_2))
    return data_conversion


# import file and set parameter
# 增加打开对应文件路径
filepath = "4.txt"
data = np.loadtxt(filepath)
data=data_conversion(data)
#增加按钮选择磁性和非磁性材料
mode = 1  # 0: nonmagnetic material, 1: magnetic material
#增加列表填写要显示的厚度
h_show = [1, 2, 3.5, 4, 5, 8,12,18]  # thickness in the RL fig, unit: mm
h = np.arange(1, 5.05, 0.05)  # thickness of the sample, unit: mm
filename = filepath.split('.')[0]

# electromagnetic parameters
frequency = data[:, 0]  # frequency, unit: Hz
e_1 = data[:, 1]  # real part of the permittivity
e_2 = data[:, 2]  # imaginary part of the permittivity
u_1 = data[:, 3]  # real part of the permeability
u_2 = data[:, 4]  # imaginary part of permeability
tan_e = e_2 / e_1  # dielectric loss tangent
tan_u = u_2 / u_1  # magnetic loss tangent

#计算要显示的RL和3维绘图RL 和阻抗匹配|Zin/Z0|
# calculate Reflection Loss(dB) and |Zin/Z0|
RL = np.zeros((len(frequency), len(h)))
Zin_Z0t = np.zeros((len(frequency), len(h)))
RL_show = np.zeros((len(frequency), len(h_show)))
Zin_Z0show = np.zeros((len(frequency), len(h_show)))
for i in range(len(h)):
    RL_temporary, Zin_Z0 = RL_calculation(data, h[i], mode)
    RL[:, i] = RL_temporary
    Zin_Z0t[:, i] = Zin_Z0
for i in range(len(h_show)):
    RL_temporary, Zin_Z0 = RL_calculation(data, h_show[i], mode)
    RL_show[:, i] = RL_temporary
    Zin_Z0show[:, i] = Zin_Z0

# plot RL
handle_1 = plt.figure(1)
handle_2 = plt.plot(frequency / 1e9, RL_show, "-", marker="s")
#选择图例
marker_sort = ["s", "o", "^", "v", "H", "<", ">", "d", "*", "h"]
for ii, line in enumerate(handle_2):
    line.set_marker(marker_sort[ii])

plt.xlim([min(frequency / 1e9), max(frequency / 1e9)])
plt.xlabel('Frequency (GHz)', fontsize=16, fontweight='bold', fontname='times new roman')
plt.ylabel('Reflection Loss (dB)', fontsize=16, fontweight='bold', fontname='times new roman')
plt.plot(frequency / 1e9, np.ones((len(Zin_Z0show), 1)) * -10, "--k")
legend_set = [f"{h_show[i]:0.1f}mm" for i in range(len(h_show))]
plt.legend(legend_set, loc='best')
plt.show()

handle_3 = plt.plot(frequency/ 1e9,Zin_Z0show)
for ii, line in enumerate(handle_3):
    line.set_marker(marker_sort[ii])

plt.xlabel('Frequency (GHz)', fontsize=16, fontweight='bold', fontname='times new roman')
plt.ylabel('|Zin/Z0|', fontsize=16, fontweight='bold', fontname='times new roman')
plt.xlim([min(frequency / 1e9), max(frequency / 1e9)])
plt.plot(frequency / 1e9, np.ones((len(Zin_Z0show), 1)), "--k")
legend_set = [f"{h_show[i]:0.1f}mm" for i in range(len(h_show))]
plt.legend(legend_set, loc='best')
plt.show()
#计算衰减常数并绘图
attenuation = attenuation_calculation(data,mode)
coefficient = u_2*(u_1**-2)*(frequency**-1)
plt.plot(frequency/ 1e9,attenuation)
plt.xlabel('Frequency (GHz)', fontsize=16, fontweight='bold', fontname='times new roman')
plt.ylabel('Attenuation', fontsize=16, fontweight='bold', fontname='times new roman')
plt.xlim([min(frequency / 1e9), max(frequency / 1e9)])
plt.show()
#计算磁响应常数并绘图
plt.plot(frequency/ 1e9,coefficient)
plt.xlabel('Frequency (GHz)', fontsize=16, fontweight='bold', fontname='times new roman')
plt.ylabel('C0/ (Hz^-1)', fontsize=16, fontweight='bold', fontname='times new roman')
plt.xlim([min(frequency / 1e9), max(frequency / 1e9)])
plt.show()
#显示特定厚度下的3维图
X, Y = np.meshgrid(frequency / 1e9, h)
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(X, Y, RL.T)
ax.set_xlabel("Frequency (GHz)")
ax.set_ylabel("Thickness (mm)")
ax.set_zlabel("Reflection Loss (dB)")
plt.show()
